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Abstract 

In this paper a new Runge-Kutta type scheme is introduced for nonlinear 
stochastic partial differential equations (SPDEs) with multiplicative trace class 
noise. The proposed scheme converges with respect to the computational effort 
with a higher order than the well-known linear implicit Euler scheme. In com- 
parison to the infinite dimensional analog of Milstein type scheme recently pro- 
posed in [Jentzen & Rockner (2012); A Milstein scheme for SPDEs, Arxiv preprint 
arXiv:1001.2751v4], our scheme is easier to implement and needs less computational 
effort due to avoiding the derivative of the diffusion function. The new scheme can 
be regarded as an infinite dimensional analog of Runge-Kutta method for finite di- 
mensional stochastic ordinary differential equations (SODEs). Numerical examples 
are reported to support the theoretical results. 
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1 Introduction 

In the last two decades, much progress has been made in developing numerical schemes 
for stochastic partial differential equations (SPDEs), see, e.g., |H El El El [HI [121 H31 HH 
[161 El [18], and an extensive list of references can be found in the review article [7]. 
In this article we are concerned with strong approximations (see Section 9.3 in [10]) to 
nonlinear SPDEs of evolutionary type. For simplicity of presentation, we concentrate 
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on the following SPDE in this introductory section and refer to Section [3] for multi- 
dimensional space case. To be precise, we consider a parabolic SPDE with multiplicative 
trace class noise as 



dX t (x) 

X t (0) - 
X (x) 



= k^X t (x) + f(x,X t (x)) 

x t (i) = o, 

■■S(x), 16(0,1). 



dt + g(x, X t (x))dW t (x), < t < T, 



Here /, g : (0, 1) xR — > K. are two appropriate smooth and regular functions with globally 
bounded derivatives, and T is a positive constant. Let H = L 2 ((0, 1), R) and let (fl, J 7 , P) 
be a probability space with a normal filtration {^t}o<t<T- Moreover, let W : [0, T] x $7 — > 
if be a standard Q- Wiener process with respect to {J~t\a<t<Ti with a trace class operator 
Q : H if. We assume that rjj,j 6 N is an orthonormal basis of if consisting of 
eigenfunctions of Q such that Qrjj = fijf]j:j 6= N. 

The problem (11. ip can be formulated in an abstract form 



dX t = {AX t + F{X t )) dt + G{X t )dW t , X = e, 



1.2) 



where A : T>(A) G H H is the Laplacian with Dirichlet boundary conditions times 
the constant k > and F : H — > H and G : H — )■ HS(Uo, H) are, respectively, given 
by (F(t>))(x) = and = f(x)) x -u(x) for all x G (0,1), v G if 

and all u G £/o- Here f/o = Q^(H) and HS(Uo, H) denotes the space of Hilbert-Schmidt 
operators from Uq to H (see Section [2] for more details). Under the assumptions above, 
the SPDE (11.2)) has a unique mild solution with continuous sample path (see, Proposition 
I2.ip . given by 



X t = S(t)£+ / S(t - s)F(X s ) ds + / S(t - s)G{X s ) dW s , P-a.s 



;i.3) 



where we denote by S(t) := e , t > the semigroup generated by the operator A. 

Now we are interested in the strong approximation problem of the SPDE (11.11) . More 
formally, we want to design a numerical approximation Y : Q — > H such that 



E 









X T -Y 




[/ 






L 7(0,1) 



\X T (x) -Y(x)\ 2 dx 



< e 



(1.4) 



holds for a given precision e > with the least possible computational effort. To imple- 
ment the numerical approximation on a computer, one has to discretize both the time 
interval [0,T] and the infinite dimensional space H. In this article we consider spectral 
Galerkin method for spatial discretization and difference method for temporal discretiza- 
tion. A simple fully discretization for ( II .ip is the linear implicit Euler scheme combined 
with spectral Galerkin method given by Y iV < M < K = p N (£) and for m = 0, 1, . . . , M - 1 

Y»>"* + hf(-,YW^+g(-,YW^xAWZ> K 



Y, 



N,M,K 
m+1 



Pi 



N 



I-hA 



-i 



(1.5) 

where we used the notations <p(-, v) : (0, 1) — > K and v x w : (0, 1) — > R given by 

((p(-,v))(x) — <f(x,v(x)), (v x w)(x) — v(x) x w(x) (1.6) 
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for all x G (0,1) and all functions v,w : (0,1) — > R and (p : (0,1) xl 4 R. Here 
and throughout this article, /i = |,M G N is the time stepsize and the increment 
AW^' K (u) := W^ +1)T (u) — W% T (u>) is given by (12.131) . The linear projection operator 

M M 

P/v : H — y H is defined by (12. lip . Here {ei}^ is an orthonormal basis of U consisting of 
eigenfunctions of A. It is worthwhile to point out that the scheme (11.51) for (jl.ip is easy 
to implement (see Figure 2 in [pj for the matlab code). 

Recently, Jentzen and Rockner [9] introduced an infinite dimensional analog of Mil- 
stein type scheme for (II. ip . given by Y ' ,K = Pn(0 and for m = 0, 1, M — 1 

= P N (s(h) (f^ K +hf(; Y^ K ) +g(., f ™) x AW^ K (1.7) 



+\{§-g) (;Y^ K ) x g (;Y^>«) x ((AW^ K ) 2 - /»J>fo 

y 3=1 



K 




Here, apart from (II. 6p . we also used the notations v 2 : (0, 1) — > R and (-^g)(-, v) : (0, 1) — > 
IR given by 

/ d \ d 

(v 2 )(x) = {v(x)) 2 , {g^9j {-,v){x) = —g{x,v(x)) 

for all x G (0,1) and all functions v,w : (0,1) — > IR. For continuously different iable 
function g : (0, 1) x R -4 R, ^g : (0, 1) x R — > R is a partial derivative of g with respect 
to the second variable. On the one hand, it is also easy to implement the scheme (II. 7p for 
the SPDE (II. ip (see Figure 3 in [pj for the matlab code). On the other hand, the scheme 
(ll.7p gives a break of complexity of the numerical approximation of nonlinear SPDEs with 
multiplicative trace class noise. For example, the scheme (11.51) for the first test example in 
Section [5] can only achieve overall convergence order |— while the scheme (II .7p possesses 
overall convergence order |— (here and below we write b— for the convergence order if 
the convergence order is higher than b — e for every arbitrarily small < e < b) . 

Note that the scheme (ll.7p can also be adapted to solve a SPDE system (II .ip with 
X t (x),£(x) G R n , and f,g : (0, 1) x R n R n for n G N. In this situation, §J in f TTTj) 
is interpreted as the Jacobian matrix of g. Therefore, to implement (JTTTJ) one needs to 
calculate the Jacobian exactly, which may be difficult, and to evaluate it at each time 
step, which may be expensive. To save computational cost in this sense, we will take a 
Runge-Kutta type scheme to avoid computing the Jacobian. Following this idea, in this 
article we aim at constructing a high strong order Runge-Kutta method for (II. ip . For 
simplicity, we only consider a scalar SPDE. But our work can be easily extended to a 
SPDE system with scalar noise as described above. 

One approach for deriving a Runge-Kutta method is to replace the partial derivative 
in the approximations (ll.7p by difference, and this leads to a derivative-free scheme, given 
by Y N ' M ' K = P N (£) and for m = 0, 1, M - 1 

Y^f K = P N (s(h) [y^ k + hf(; Y^ K ) + g (, y W) x AW^ K (1.8) 
g(;Y^ K + Vhg(.,Y^ K )) - g(- ,Y»> M > K )} x ((AW^ K ) 2 - h f>(^) 2 



j=i 
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where the function g ( •, v + \/hg{-, v) ) — g(-,v) : (0, 1) — > R is defined by 




for all functions v : (0, 1) — > R, g : (0, 1) x R — > R. A natural question thus arises as to 
whether such replacement maintains the high convergence order of (11. 7p . In this paper 
we give a positive answer and prove that the new scheme not only maintains the high 
convergence order, but also reduces computational cost. Similarly to the scheme (jl.7p . 
the numerical method f ll.Sp can be simulated quite easily (see Figure [T] in Section \5\ for 
the implementation code). 

Now we take a closer look at schemes fll.7j) and (II. 8p . For each step, the Milstein type 
scheme fll.7j) requires one evaluation of /, one evaluation of g and one evaluation of the 
partial derivative -j^g. In contrast, the new Runge-Kutta type scheme (11.81) needs one 
evaluation of / and two evaluations of g, but no evaluation of the partial derivative -^-g 
at each step. Thus the Runge-Kutta type scheme (ll.8p is easier to implement than the 
Milstein type scheme (II. 7p . The main result (Theorem 12 .ip shows that the new scheme 
(ll.8p maintains the high convergence order of scheme (11.71) . Numerical results in Section 
[5] demonstrate that the schemes fll.7[) and (I1.8P produce nearly the same approximation 
errors. Even if we neglect the effort for the calculation of the partial derivative, the 
runtime for one path simulation of the Runge-Kutta scheme (II. 8p applied to (ll.ip is less 
than that for simulating the scheme (11.71) . Take the first test problem in Section [5] for 
example, in the case of iV = 256, one path simulation of the Runge-Kutta scheme (ll.8p 
costs 60.312000 seconds, while it needs 77.016000 seconds to simulate the scheme (11.71) (see 
Tabled]). This occurs due to the fact that evaluation of the partial derivative (^g)(x,y) 
costs more time than evaluation of the function g(x, y). Summarizing, the Runge-Kutta 
type scheme (ll.8p is easier to implement and needs less computational effort than the 
Milstein type scheme (ll.7p . Moreover, the derivative- free scheme (ll.8p can be regarded as 
an infinite dimensional analog of the Runge-Kutta method (6) in [1] for finite SODEs. To 
the best of our knowledge, this is the very first paper to introduce a Runge-Kutta method 
for nonlinear SPDEs with multiplicative noise. We also mention that constructing higher 
order Runge-Kutta methods and developing more systematic way to derive Runge-Kutta 
methods for SPDEs are of great interest and will be our future work. 

The rest of this paper is organized as follows. In the next section, we put everything 
into an abstract framework and state the main convergence result of this article. In 
Section [3] we give examples fulfilling the assumptions in the previous section. A detailed 
proof of the main convergence result is elaborated in Section HI Finally, we illustrate how 
to implement the proposed scheme and present some numerical examples to support our 
theoretical results. 

2 Abstract framework and main result 

In this section we focus on the abstract framework (II. 2p and adopt the following setting 
and assumptions. 

Let (H, (•, -) H , \\-\\ H ) and (U, (•, ') v , \\'\\ v ) be two separable Hilbert spaces. By L(U, H) 
and L( 2 \U, H) we denote the space of linear bounded operator from U to H and from 
U xU to H, respectively. For short, we write L(H) instead of L(H, H). We also introduce 
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a space of Hilbert-Schmidt operators. An operator T G L(U, H) belongs to the Hilbert- 
Schmidt operator space HS(U, H), if for any orthonormal basis {ipk}kLi of U the sum 

oo 

l|r|lHS(E/,io := H^^lln 

fc=l 

is finite and is independent of the choice of the orthonormal basis. The quantity ||r||/fs(r/,H) 
is called the Hilbert-Schmidt norm of V. Similarly, we can define an Hilbert-Schmidt op- 
erator space HS^(U, H) from U xU to H . We refer to Chow [2], Da Prato and Zabczyk 
[3], Prevot and Rockner [15] for details on these spaces and their properties. 
Moreover, let Q G L(U) be symmetric, nonnegative and with finite trace, i.e., 

Tr(Q) < oo. (2.1) 

Suppose that J is a finite or countable set, and let {i]j)j<=j C (7 be an orthonormal 
basis of U consisting of eigenf unctions of Q : U —*U such that Qr/j = G J. 

Denote by (U , (-, -) Uo , ||-||;y ) the separable Hilbert space U := Q^{U) with (v,w) v = 

(Q~^v,Q~^w) v for all v , w G Uq, where is the pseudo inverse of (see, e.g., 
Section 2.3.2 in [15]). We can obtain that 



\T\\hs(u ,h) 



ToQi 



for TeHS(U ,H). (2.2) 

HS(U,H) 



Further, we assume that (f^J 7 , P) is a probability space with a normal filtration 
{^~*}o<i<T and W : [0, T] x Q, — > U is a standard Q- Wiener process with respect to 
{Ft\o<t<T and has the representation [151 Proposition 2.1.10] 

W t {u) := £ y/fiPtMnj, (2.3) 

where (Pi)jeJ,nj^o f° r £ G [0>^1 are independent real-valued Brownian motions on the 
probability space (fl, J 7 , {J r t}o<t<T ) P)- Now we make the following assumptions. 

Assumption 2.1 (Linear operator A). LetX be a finite or countable set and let 

be a family of real numbers with infj e j \ G (0, oo). Further, let {eA^x be an orthonormal 

basis of H and let —A : T>(—A) C H — )■ H be a linear operator such that 

- Av = J2\ i (e i ,v) H e i (2.4) 

iei 

for allveV(-A) := {v G H\ £ ieI \X t \ 2 \ (e u v) H \ 2 < oo} . 

Here and below we denote by S(t) := e At ,t > the semigroup generated by the 
operator A. By V r := V ((— A) r ) , r > equipped with the norm ||f || v := ||(— A) r v\\ H we 
denote the M-Hilbert spaces of domains of fractional powers of the linear operator —A. 

Assumption 2.2 (Drift coefficient F). For (5 G [0,1), we assume that F : Vg — )■ H is 
a twice continuously Frechet dijferentiable mapping with sup veV \\F'(v)\\ L ( H j < oo and 
with sup„ ev> \\F"{v)\\ L(2){v ^ H) < oo. 
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Assumption 2.3 (Diffusion coefficient G). Let G : Vp — > HS(Uq, H) be a twice contin- 
uously Frechet differentiable mapping such that sup veVfj \\G'(v)\\nH,HS(u ,H))< oo and 
snp veVfj \\G"(v)\\ L( 2) {v ^ HSiUo!H)) < oo. Moreover, let a, c E (0,oo), 6,-& E (0, |) with 
j3 < 5 + \, and 7 E [max (5, 5 + ~). Suppose that G(V$) C HS(U , Vs) and 

ii^Wo^^^i+nu, (2.5) 

G(v) - G\w) G(w) \\ HSi2)(U0)H) <c\\v-w\\ H , (2.6) 
-A)-°G(v)Q- a _ ^ < c (l + |M|„) (2.7) 



HS(Uo,H) 

hold for all u E Vs and v,w E V 1 . Furthermore, let the bilinear Hilbert- Schmidt operator 
G'(v)G(v) E HSW(U , H) be symmetric for all v E Vp. 

Note that the operator G'{v)G{v) : Uq x U — > H given by 

(G'(v)G(v)) (u, u) = (G'(v) (G(v)u)) (fi) (2.8) 

for all u, u E Uq is a bilinear Hilbert-Schmidt operator in HS^ 2 \Uq, H) for all v E Vp. 
The assumed symmetry of G'(v)G(v) E HS^ 2 \Uq, H) thus reads as [9j Remark 1] 

G'(v) (G(v)u)) (it) = (g» (G(u)u)) (u) (2.9) 

for all u,u E Uq and all t> G Vg. We also mention that (12. 9p is the abstract (infinite 
dimensional) analog of the commutativity condition (10.3.13) in [10]. Although the com- 
mutativity condition (10.3.13) in [10] is seldom fulfilled for finite dimensional SODEs, 
(123)]) is naturally met for SPDE (fTTTl) . 

Assumption 2.4 (Initial value £). Let £ : fl — >• 6e an Tq/B (Vy) -measurable mapping 
wUhE\\£\\y < oo. 

We remark that a mapping £ : fl — > V y is called J-q/B{V^) measurable if it is a mea- 
surable mapping from the measurable space (fl, J-q) to the measurable space (Vy, B(Vy)). 
Here B(V y ) denotes the Borel a-field of V 1 . The assumptions above are sufficient to 
guarantee the existence of a unique mild solution of the SPDE (II. 2p [8; Theorem 1]. 

Proposition 2.1 (Existence of the mild solution). Let Assumptions \2. l\2.J\ and condition 
(12. ip be fulfilled. Then there exists an up to modifications unique predictable stochastic 
process X : [0, T]xfl — > V 1} which fulfills sup E||X t || v < oo ; sup E || Gr(X t )|| HS ^ Uo y^ < 



oo and 



X t = S{t)i+ [ S(t-s)F(X s )ds+ [ S(t- s)G(X s )dW S} F-a.s. (2.10) 
Jo Jo 

for allt E [0,T]. 

Let (X A r) 7VgN and {JkIkg* be sequences of finite subsets of I and J , respectively. For 
N EN we define the linear projection operators Pjv : H — > H, by 

Pn{v) := X)< e *' v >H^ veH - ( 2 - n ) 
6 



Furthermore, for all K G N we define Wiener processes W : [0,T] x Q — >■ Uq by 

Wf (w) := £ v % i G [0,T], w G fi. (2.12) 

Here {^)j€J,n-^Q are independent real- valued Brownian motions on the probability space 
(fi, J 7 , {J r t}o<t<T, IP)- We also use the notations Aft^u) := /3f (tu) — Pf m (u) and 

: = WfcfOT (g;) - ( W ) = J2 sfaWMVi (2-13) 

for all uj G f2 and m = 0, 1, . . . , M — 1. 

Subsequently, we formulate the schemes in the introduction part in abstract form. 
Then the scheme ( ll.7p for the problem ( 11. ip can be formulated in an abstract scheme for 
the abstract problem (II. 2p . given by Y ' ' = Pn{£,) and for m = 0, 1, . . . , M — 1 

yN^K = p N ^ {h) (yN,M,K + fc jr (y w) + (?(yW)A^ (2.14) 

+ \<? (g(y^)aw^)aw^ - \^ G '{^ M ' K ) 



\M,K 



Vj Vj 



Assume that Assumptions 12.1112.41 are all satisfied, Jentzen and Rockner [9] have es- 
tablished the strong convergence result for the scheme (12.141) (see Theorem 1 in [9j). In 
later development, we show that this convergence result also holds for the Runge-Kutta 
type method proposed here under some additional conditions. First of all, we formulate 
the new numerical scheme (ll.8p in an abstract form, given by Y Q N,M ' K = Pn(0 and for 
m = 0,l,...,M-l 



Y *M* =Pn ( S (h) \YW>* + h F(Y r ^' M ' K ) + G (Y^ M ' K ) AW^ (2.15) 
+ \gG(Y^ k , h) (AW^ K , AW^ K ) ^GG{Y^ k , h) ( Vj , Vj ) 

where GG(v, h) : Uq x Uq — > H is a derivative-free bilinear operator to approximate the 
bilinear operator G'(v)G(v) : Uq x Uq — >■ H in ( I2.14p . We define the remainder bilinear 
operators GG\(v, h) : Uq x Uq — > H by 

(gGx(v, hj\ (u, u) := (GG(v, hj\ (u, u) - (g'{v){G{v)u)} (u). (2.16) 

Note that Assumptions I2.1H2.4I come from [9], which are used to ensure the high 
strong convergence order of the Milstein-type scheme (11.71) . In the next section, concrete 
conditions are given for concrete parabolic SPDE to promise these assumptions. To guar- 
antee the high convergence order of the Runge-Kutta scheme, we shall impose additional 
conditions on the two bilinear operators GG(v, h) and GG\(v, h) as follows. 
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Assumption 2.5 (Approximation operator). Suppose that for any v,w E.Vp there exists 
a constant Co independent of h such that 

\\GG(v, h) - GG(w, h)\\ 2 HSi2)(Uo H) <Q\\v - w\\ 2 H , (2.17) 
||C?Gi(t;,/i)||^ (l7o|jB) <Co/i(l+||t;||y. (2.18) 

In the next section we will validate the imposed conditions (I2.17p -f l2.18l) for parabolic 
SPDEs (see Proposition 13.11) . Armed with these assumptions, we are now to give the 
main result of this article. 

Theorem 2.1 (Main result). Suppose that Assumptions Iff. llTOl and ( 12. ip are fulfilled. 
Then there is a constant C > independent of h such that 

sup (E\\X mh -Y»> M >%f <c(( inf A,r 7 +( sup ^)\m~^-^A. 
o<m<M v ' \ ^ iex\z N f V jeJ\J K J J 

(2.19) 

The detailed proof is postponed to Section HI Theorem 12.11 indicates that the approx- 
imation error in (12.191) is composed of three parts. The first term (inf i^x\x N Aj) 7 arises 
due to spatial discretization. The second term (^Pj & j\j K comes from truncation 
of the expansion of the noise W t . The third term m~ uiM2 ^~^^ corresponds to the 
temporal discretization error. 



3 Parabolic SPDEs 

In this section we give concrete parabolic SPDE examples falling into the abstract frame- 
work in Section [2J Let d G {1,2,3} and let H = U = L 2 ((0, l) d ,lR) be the Hilbert space 
with the scalar product and the norm, respectively, given by 



(v,w)h = / v(x) x w(x)dx and ||t> \\h = I / \ v ( x )\ dx 
J(o,i) d \J(o,i) d 

For the continuous function v : (0, l) d — > K, we define two norms 

\\v ||c((o,i) d « := SU P \ v ( x )\ 

26(0, l) d 

and 

IMIo((o,i) d « := su p \ v ( x )\ + su p I, _ .i r — , 

a;G(0,l) d x,y£(0,l) d ,x^y \\ x VuRd 

where the Euclidean norm \\x\\^d := (|xi| 2 + ... + |xd| 2 )5 was used for x = (xi, Xd) G K d - 
First of all, we assume that for some constants < p < l,c > 0, the eigenfunctions 
Vjjj G J7" of the covariance operator Q are continuous and satisfy 



sup \\Vj\\c((o,i) d m < c, ^2^j\\Vj\\ 2 cP( { o,iy ;R) < c (3.1) 
Here and below c is a generic constant, which may be different in different places. 
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For the linear operator A in Assumption I2.1[ let I = N d and A := kA with k > 
be the Laplacian times a constant with Dirichlet boundary condition, i.e., Av = kAv = 
k \^2j=i Ji?) v ^ or v ^ D{—A). Then (12. 4p in Assumption 12. II holds with 

d d 

2a Y\ sin(ij7TXj), A, = kir 2 ^""^(i?) 2 



for all x = (xi, x d ) G (0, l) d and all i = z^) G N d . Here we set = {1, N} d . 

For the drift coefficient F in Assumption 12.21 set 

/? = -, for d=l,2,3, 
5 

and let / : (0, l) d x 1 -> 1 be a twice continuously differentiable function such that 



(o,i) d 



\f(x,0)\ 2 dx <c, 



d n 
dy n 



fifay) 



< c, n— 1,2 



for x G (0, l) d , y G R. Then we define the operator F : Vg — )■ if as 

(F(v))(x):=f(x,v(x)) 

for all x G (0, l) d and all v G Vg. 



(3.2) 



(3.3) 



For the diffusion coefficient G in Assumption I2.3[ let g : (0, 1) x R — >• R be a twice 
continuously differentiable function with 



< c, 



d" 



dy r 



9) fay) 



< c, 



9) fay) 



and 



d_ 

dy l 



g)(x,y)g{x,y) 



d_ 

dy l 



d_ 
dx' 



g)(x,z)g{x,z) 



< c, n = 1,2, 



< c\y — z\ 



(3.4) 



(3.5) 



for x G (0, l) d , y,z G R. Here || • Hl^r) is the usual operator norm. Then let the operator 
G : Vp ->■ #£([/o, F) be given by 

(G(f )m)(x) := g(x, f (x)) x -u(x) 

for all x G (0, l) d , u G Vp and w G t/o C U = H . Therefore the bilinear operator G'(v)G(v) 
in the abstract scheme (12.14p is here given by 

G'(v)G(v)j(u,u)(x) = ^—gj(x,v(x))g(x,v(x)) x ufa) x ufa). (3.6) 

In regard to the initial value in Assumption I2.4[ let xo : [0, l] d — > R be a twice 
continuously differentiable function with Xo|a(o,i)<* = 0. Then let the initial value be given 
by £(u>) = Xo for all u G Vt. 

It is shown in [9j Section 4] that the linear operator A, the drift coefficient F, the 
diffusion coefficient G and the initial value £ defined as above satisfy all conditions in 
Assumption I2.1H2.4[ except for (12. 5p and (12. 7p in Assumption 12. 3\ which will be verified 
for some concrete examples in Section [5j 
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In the setting above, SPDE (II. 2p reduces to a parabolic SPDE as 



dXt(x) 



k[J2-^)x t (x) + f(x,X t (x)) 



dt + g(x, X t (x))dW t {x) (3.7) 



with X t \ 9 r 0jl \d = and X (x) = xo(x) for x G (0, l) d and t G [0, T]. For ( 13. 7p . linear 
implicit Euler method and Milstein type method take the same form as (II. 5p and (II. 7p . 
respectively. If we introduce a bilinear operator GG(v, h) : Uq x Uq — >■ i/ approximating 
the bilinear operator G'(v)G(v ) in the scheme (I2.15p . given by 



GG(v, h) ) (w, = — y= \ g [x, v(x) + Vhg(x, v (x)) ) — g(x, ) x x u(x) 



(3.8) 



for all f G V/3, then the scheme (I2.15P reduces to the concrete scheme (II. 8p . 

Apart from (I2.5P and (12.71) . one also needs to verify Assumption 12.51 for (13. 7p . 



Proposition 3.1. Suppose that the bilinear operators G'{v)G{v) and GG(v, h) are given 
by ( 13. 6 p and (13. 8p . respectively. Then the operator GG(v,h) and the remainder opera- 
tor GG\(v,h) given by (12.161) fulfill the conditions in Assumption \2.5[ provided that the 
conditions (13. ip and (13. 4p hold. 

Proof. From (I3.8P we have 



\\GG{v, h) - GG(w, h)\\ 2 Hsm(UQtH) = Mi GG{v, h)^,^) - GG(w, h)^,^) 

i,ieJ 



-y 
h ^ 

i,iej 



(o,i) d 
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x, v{x) + Vhg(x, v(x)) ) — g(x, v(x)) J x r]i{x) x r)j(x) 



x, w(x) + Vhg(x, w(x)) ) — g(x, w(x)) J x r)i(x) x T]j(x) 



dx 



<-y 

i,iej 



x, v{x) + Vh ■ g(x, v(x)) ) — g[x, w(x) + Vhg(x, w(x)) 



+ 



(o,i) d 



g(x,v(x)) - g(x,w(x)) 



dx X |hj|lc((0,i) d ;lR) X 



2 

C((0,l) a 



(3.9) 



Using ( 13 .4p shows that 



\g(x, v(x)) — g(x, w(x))\ < c \v(x) — w(x)\ 



(3.10) 



and that 



x, v (x) + Vhg(x, v(x)) J — g( x, w(x) + Vhg(x, w(x)) 



<c A 



v(x) + Vhg(x, v(x)) — w{x) — Vhg(x, w(x)) 



<2c 2 \v(x) -w(x)\ 2 + 2c 2 h 



g(x,v(x)) - g(x,w(x)) 



<2c 2 (l + Tc 2 )\v(x) -w{x)\ 2 , 



(3.11) 
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where we also used the fact that h < T. Inserting (13.101) and (13.111) into (13.91) yields 

\\GG(v,h)-GG(w,h)\\ 2 HSi2){Uo>H) 

<t Yl I c 2 (3 + 2Tc 2 ) \v(x) - w{x)\ 2 dx x II^H^i^r) x IMIc((o,i)<* ; i 

11 i,jej V J (°^ d 



< 



2c 2 (3 + 2Tc 2 
h 



W ( 



i,jej 



SUP ||77i|| C((0 ,l)rf;I 



\v — w\ 



H 



2c 6 {3 + 2Tc 2 ){TrQy 
h 



12 
\H- 



Now the estimate (12.171) in Assumption 12.51 is validated on choosing C > 2c 6 (3 + 
2Tc 2 )(TrQ) 2 . 

For the second estimate (12.181) . we use Taylor's formulae in (13.81) and derive that for 
all v G Vp the remainder operator GGi(v , h) : Uq x Uq — y H given by fl 2 . 1 6 [) satisfies 

GG 1 (v,h))(u,u)(x) = (GG(v,h))(u,u) - (g'(v)(G(v)u)){u) 



1 (d 2 



79 



x, v(x) + rVhg(x, v(x)) ) (1 — r) g 2 (x, v{x))dr x u{x) x u{x). (3.12) 



'o -dy 2 

Combining (I3.4[) and ( 13 . 1 2[) and taking (12 .2p into account show 
\\GGi{v,h)\\ 2 HS(2){Uo . H) = Y HiHj\\GGi( v i h )(Vi,Vj)\\ 2 H 



h \ 



1 , Q2 



dy 



;9 



x, v(x) + rVhg(x, v(x)) 



x (1 — r) g 2 (x, v (x))dr x rji(x) x rjj(x) 



<c 2 h V ^ / g{ 



dx 



(x, v(x)) 



dx x 



* IIC((0,l) d ;R) X 11 !' ,; 1 >•'■ 



2 

C((0,l) d 



Using the elementary inequality (|a| + \b\) p < 2 P 1 (\a\ p + \b\ p ) for p > l,a, b G R, and 
rtTOjt gives 

|^(x,t;(x))| 4 < (b(a;,«(a;)) - #(x,0)| + |#(:r,0)|) 4 < (c|t;(x)| + c)* < 8c 4 (|t;(x)| 4 + 1 
Thus for all v G Vg, 



||GGi(u,/0||L? (a)(Ob . H) ^ 8c6/l S ^ ( sup 



<8c 10 (TrQ) 2 M|| W || 4 yfl + l), 



4 



4 

C((0,l) d ; 



'lll,4((0 ; 



where (13.11) and the fact were used that Vp C L 4 ((0,l) d ;M) continuously for p = jj by 
Sobolev embedding theorem. The proof of Proposition 13 . 1 1 is complete. □ 

So far, all conditions in Assumption I2.1H2.5I have been verified except conditions (12. 5 j) 
and (12.71) in Assumption 12.31 We mention that condition (12. 5p originally comes from [8], 



11 



where (12. 5p is a key ingredient to promise higher spatial and temporal regularity of mild 
solution. The condition ( 12 .7p is needed to estimate approximation error due to truncation 
of the expansion of the noise Wt- In the last section, we will obtain these two conditions 
for two concrete examples. 



4 Proof of Theorem 12.1 



First of all, we rewrite the numerical solution ( I2.15P in the following form 

An— 1 

= S(mh)P N (0 + P N P£ S((m - l)h) F (y™*) ds 

E / (4.1) 



1 / to— 1 
\ Z=0 



m— 1 



1=0 jejx 
Likewise, the numerical solution (I2.14p satisfies 

'm-\ r (l+l)h 



(I I C — -1- n ( J. I fl 

Ej lh s({m-t)h)F(lr»> u >«)d8 

£ / ^((m-O^G^^die (4.2) 

;_n J Ih I 



1=0 
f m—\ 



+ \ Pn ( E 5 (( m " 0^)G , (f^ M ' K )(G'(^ M ^)A^ M ^) (aw { - 

V 2=0 

, / TO— 1 



Moreover, the exact mild solution of the SPDE (11.21) can be rewritten as 

rmh pmh 

X mh = S(mh)£ + / S{mh - s)F(X s ) ds + S(mh - s)G{X s ) dW s (4.3) 
Jo Jo 

m ~ l Ml+l)h m - 1 r(l+l)h 

= 5(m/i)^ + E / S(mh- s)F(X s )ds + J2 / S(mh - s)G(X s ) dW s , 
1=0 -> lh i=o 

and thus 



E / S(mh - s)F(X 8 ) c 
i=o Jlh 

/to-1 p{l+l)h \ 

+ Pn [J2J i S(mh-s)G(X s )dW s j. (4.4) 
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To estimate E \ \X mh — Ym' M,K \\w we need two auxiliary processes Z^' M,K and 
Define Z% M > K by 

tm—l p(l-\-X)h ^ 
J2j ih s((m-l)h)F(X lh )ds 

J2j ih s((m-l)h)G(X lh )dW s K \ (4.5) 



1=0 

+ 2 Pn (Y, S {( m ~ l)h)GG(X lh , h){AW^ K , AW^ K ^ 
V 1=0 

, / m— 1 

Note that (14.1 p coincides with (14. 5 p with Y^> M ' X replaced by X^. Similarly, replacing 
yN,M,K j n (|4 2 p by we introduce the process Z^' M,K , given by 

(m— 1 ,>(7+i)/i ^ 
S((m-I)fc)f(4)* 

J2j ih s({m-l)h)G{X lh )dW«\ (4.6) 



+ -Piv J] 5 (( m " l ) h ) G '(Xih) (G{X lh )AW^ K ) [AW l 
\i=o 

fm—l 
^S((m-l)h) Y,^G'{X lh )(G{X lh )r ]j 



M,K 

I 



h 

Pv( 

«=0 ' ' j£j K 

H ^0 

Armed with these notations, now we start the proof. Employing the elementary 
inequality (ai + a2 + a^) 2 < 3(|ai| 2 + |a 2 | 2 + 1 0,3 1 2 ) , a\, a 2 , a 3 G M shows for m — 0,1, M 

E\\X mh -Y»> M > K \\ 2 H (4.7) 

I 2 >11W\\P.,(V .\_7N,M,K\\Z ,o W \\7N,M,K_ v N,M,K\\2 



\H ' 



<3E \\X mh - P N (X mh )\\ z H + 3E \\P N (X mh ) - Z% M > K \\ H + 3E \\Z% M > K - ^ 

Below we will estimate the three terms in (14. 7p step by step. First let R > be a real 
constant such that 

\\F'(v)\\ H h) < R, \\G'(v)\\ L{H>HS{Uo>H)) < R, 
E\\(-AyX t \\ 2 H = E\\X t \\^<R, E\\X t \\^<R 

for all v G Vp and all t G [0, T]. Due to Assumptions 12 .1112 .41 in Section [2] and Proposi- 
tion 12. H such a real constant exists. 

For the spatial discretization error E \\X m h — P/v(X m / l )|| 2 ? , using (14. 8p we derive 



E || X mh — P N (X mh ) \\ 2 H — E || (/ — P N ) X mh \\ 2 H 



<M-A)-^I-P N )\\*xE\\X.. 



L(H) " " W^rahWy^ 

«GI\Ijv 



-2-y 

<R\\(-A)-^(I-P N )\\ 2 L{H) = R( inf Aj ) (1.9) 
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To estimate E\\P N (X mh ) -Z%> M ' K \\ 2 H , we need the estimate E\\P N (X mh ) - Z^' M > K \\ 2 H . 

Lemma 4.1. Under Assumptions \2. 1\ 2.4 and (12. ip . there exists a constant Ci, indepen- 
dent of h, such that 



sup E 

0<m<Af 



p N (x mh ) - z»> M > K 2 <d f ( sup + M-^T-nnA . (4.io) 



Proof. To establish the convergence result for the Milstein type scheme fl2.14[) . (I4.10p 
has been obtained in [9] (see Section 5 in [9] for the details). □ 

Lemma 4.2. Suppose that Assumptions \2.1W2.5\ and the condition (12. ip are fulfilled. 
Then there exists a constant C 2 , independent of h, such that 



sup E 

0<m<Af 



Z 



N,M,K 



2 / \ 2a 

< C 2 SUp [L j + 
\j£j\J K J 



H 



Co 



^fmin(4( 7 -/3),2 7 ) 



(4.11) 



Proof. Using ( I2.16p . we derive from ( 14. 5 j) and (14 .6p that 

_ / m— 1 
\ 1=0 

, / m—l 

~2 Pn [ zZ S {( m - l ) h ) zZ ^ GG ^ X ^ h )^^) 



(4.12) 



1=0 



H ^o 



Therefore 



td (y \ yN,M,K 



-Pjsr(X m h) — Z 

f m—l 



N.M.K 



1 / !" . 

- -Pn [Y, S (( m ~ l)h)GG 1 {X m h)(AW t M ' K , AW t M ' K ) 

\ 1=0 

+ 2 Pn [ z2 S (( m - l ) h ) z2^ GG ^ Xlh ' h ^,Vj) 



and thus using the elementary inequality (at + a 2 ) 2 < 2(|ai| 2 + | a-2 1 2 ) , a-i, a 2 G E yields 



E 



)-z 



N,M,K 
m 



<2E 



H 



Pn {X m h 

)-Z 



N,M,K 
m 



+ -E 
2 



m—l 



J>((m - l)h)GG x (X lh , h) {AW^ K , AW, 

1=0 

... i 

hJ2s((m-l)h) HGGiiX^h)^,^) 



m—l 



:=2E 



H ^o 

2 

p (y \ yN,M,K 

H 



(4.13) 
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where the fact that ||P/\HI# < \\ v \\h was also used. Due to (I4.10p in Lemma I4.1[ it 
remains to estimate J\. Inserting the representation (I2.13P and using bilinearity of the 
operator GG\ give 



1 TO— 1 

h =-e|| J2 s (( m ~ l ) h ) ( GG i( x ih, hj) ( E VWAPfa, E v^ A A 



1=0 

m—l 



jeJ K 



-/ij]5((m-l)/i) E y/jr i y/Jj~5i j ^GGi( y Xi h ,h)^( y r] i ,r] j ) 



-E 



2=0 
m—l 



E E v^V^^m - 0^) /i)) fa, ^) x (A#A# - ^/i) 



H 

(4.14) 



where <5y = 1 for z = j and 5jj = for i j. For simplicity of notation, we denote 

xI J = y/J^iyJWjS (im - l)h) (GGi(X lh , h)\ fa, tj,). 
Then we can rewrite (I4.14p as 



m—l 



m—l 



ii=0 h,ji&JK 
Mil .Mji 7^0 



^2=0 l2,j2&JK 



ii 



m—l m—l 



2 E E E E E ( (xSi* i^TMi - s ^ h ) wM: - s ^ h ) 



h=0 h=0 iij'ie Jk i2,j2Gjfc 
Mil 'Mil ^° Mi 2 >Mj 2 ^° 



In the case that h h, without loss of generality we set l\ < 1%. Using the fact that 
(Af3l 2 ) i( zj K ^ i7 L are independent of Tt x D , E(A/3/) = and xf £ -^tj shows that 



E 
=0. 



Here the last step follows by the obvious fact that E (A/3/ 2 A/3/ 2 2 — 5 i2 j 2 ti) = 0. In the 
case that l\ — l 2 — I but that ^ {12,32}- Using the mutual independence of 

A/3/, A/3/, i 7^ j and the fact that (Af3 t l ) ie j K ^ iJ i are independent of J- tl , E(A/3/) = and 
Xi' 3 £ -T 7 *! shows that 



E 



XT, XZ' 32 > TT mi^ll - Knh) (A/3- A/3/ 2 - 5 hh h) 



H 



^{x\ 1 ' 3 \xT' n ) H HW^l 1 -5 n nh) (A/3- A/3f - 5 J2 , 2 / i ) 
=0. 
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Here the last step follows since 

E ((A/?; 1 A/3/ 1 - 5 hjl h) (A/3; 2 A/3/ 2 - 8 i2h h)) = 

in the case that {i±,ji} 7^ {12,32}- Hence, using the fact that x]' 3 £ J~t x and that 
(A/3;)j e j x>Ati ^o are independent of we have 

1 m—l 



=0 ijejjf 



m—l 



<EE e (A/3; a/3/) + /> 2 E E £ M" 



=0 ijejif 



m—l 



m—l 



E E [* 

z=o ijeJk 



S((m-l)h)( GG x (X lh , h) ) (^T iVl} ^JT ]V] ) ^ E| A#A# 



i A «3|2 



m—l 



+ /i 2 ^E ^ js^(m-l)h^(GG 1 (X lh ,h)y^r ]i , y /JI]r 1jj 



1=0 \ ijeJk 



H 



Due to the fact that HS^Hk < > for all v G -ff and that E|A/3;A/3/| 2 < 3/i 2 

for all i,j G J7ft-, and using (I2.2p . (12. 18[) in Assumption 12.51 and (14 .8p we derive that 

2 



m—l 



J!<3/i 2 ^E ]T 



1=0 \ i,j€j K 
m—l 



GG 1 (X lh ,h))(y/JT i r] i , y/i^rjj) 



+ h 2 J2 E \\ GG i( X ih,h) 
1=0 

m—l 

KAh^AGG^X^h) 



HS( 2 ) (U ,H) 



1=0 
m—l 



HS( 2 )(U ,H) 



<4C /i 3 ^(l + E\\X; 



4 v 



1=0 



< 



AC Q T 3 {1 + R) 
M 2 



Plugging (I4.10p and the preceding estimate into (14.131) gives 



E 



-z, 



N,M,K 
m 



2 / \ 2Q 2C 1 + 4CoT*(1 + R) 

< 2Ci sup fij + 

# \jeJ\Jx J 



(4.15) 



(4.16) 



j^min(4( 7 -/3),27) 

Consequently, (14. lip is derived on choosing C2 = 2Ci + 4cT 3 (l + R). □ 

Lemma 4.3. Suppose that all conditions in Assumptions \2.1\\2.5\ are fulfilled. Then there 



exists a constant C3, independent of h, such that 

m—l 

jg. \\gN,M,K _ Y N ' M ' K \\ 2 

1=0 

holds for all m = 0,1, ... , M. 



N,M,K 



(4.17) 
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Proof. Using the elementary inequality (ai+a 2 +a 3 ) 2 < 3(|ai| 2 +|a2| 2 +|a3| 2 ), a±, a 2 , 03 G 
gives 



U7 II n,N,M,K V N,M,K\\ 2 
^ H^m ~~ r m ||// 



<3E 



m-l ,(H-l)h 

W 5((m - l)h) (F(X lh ) - FiY^)) ds 
1=0 Jlh 

™~\ r(l+l)h . 
£ / 5((m-/)/,)(G(X z/l )-G(^ M ' K )) 
z=o 

m—l 

J>((m- Z)/i) (GG(Jf, h , /i) - GGIyW'*, fc)) (AW t M ' K , AW t M ' K ) 
1=0 

1 



+ 3E 



+ 7 E 
4 



H 



■ — J 2 + J3 + Jii 

where we also used the fact that \\PnV \\h < Ik 1 1 if- For J 2 , one can derive that 



(4.18) 



' m—l 



J 2 <3Mh 2 £ E |K( m - 0^) (^(^ifc) " ^(V' 
V 1=0 

(m—l 
J2^\\F(X lh )-F(Y^ 
1=0 

m—l 

<3TR 2 h^2E\\x lh -Y t N ' 



\M,K 



\M,K 



1=0 



(4.19) 



where (14. 8 p and the fact that 
that 



L(H) 



< l,t > were used. For J 3 , one can derive 



m-l r (l+l)h 

J 3 <3 W E\\s((m-l)h) {G{X lh )-G{Y l NM - K )) 
1=0 



HS(U ,H) 



ds 



m—l 



<3R 2 h^E llXth-Y*' 



M,K 



1=0 



(4.20) 



where the isometry property, (14.81) and the fact that ||5'(f)||/ J (//) < l,t > were again 
used. 

Now it remains to estimate J4. In a similar way as estimating J\ and using the 
condition ( 12 .17p in Assumption 12.51 and (14. 8p . one can obtain that 



m—l 



J 4 <6h 2 £ n\GG(X lh , h) - GG(Y l N ' M,K , h)\\% s(2){Uo 



H) 



1=0 
m—l 



<&c hY,nxih-Y l 



N,M,Kn2 

\\h- 



(4.21) 



1=0 
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Now, inserting (ETT9]) . K2Q\i and (OT]l into ( gJED gives the desired estimate (Q7|) . □ 
Now we return to (14.71) . With the estimates (I4.9p . (14.111) and (14.171) at hand, we derive 
from < K7} that 



m— 1 



\M,K 



< \\X nill Y^ K \\ 2 H <3C 3 h J2 E \\ X ^ - Yi N ' f 

-27 / x 2,i 



Z=0 



ii 



(4.22) 



+ 3i?( inf A, 

i£l\I N 



3C 2 [ sup fij + 



j&J\J K 



3Co 



^min(4( 7 -/3),2 7 ) 



Finally, Gronwall's lemma gives the main result (I2.19p and the proof of Theorem 12.11 is 
complete. 



5 Numerical experiments 

In this section we will first illustrate how to implement the scheme introduced in this 
work and then present some numerical results to support our theoretical assertions. 



5.1 Implementation 



For simplicity of notations, here we only consider the new scheme (ll.8p for one dimen- 
sional space case (II .ip and one can adapt the following implementation to handle multi- 
dimensional space case and other schemes. Using the notation Cm : (0, 1) — > R given 
by 



Ux) =Y^ K {x) + h f(x, Y^ K {x)) + g (x, Y^ K {x)) x AW™> K (x) 

9 ( x, Y^ K {x) + Vhg(x, Y^ K {x))) - g(x, Y^ K {x)) 



+ 



2y/h L 



A' 



(5.1) 



(Aw^ K (x)r-hj2^(v^)f 

for m = 0,1, ...M — 1, x 6 (0, 1), the scheme (II .8p can be rewritten as 

N N 

YZK* = Pn (S(h) Cm) = £ (e Ah Cm, e 3 ) H e 3 = £ e~ x > h (C m , e,) H e, 



(5.2) 



3=1 



for m = 0, 1, M — 1 and 



N 

Y W = P N {i) = Y.^)nZ3- 
3=1 

Here H = L 2 ((0, 1); M), A is the Laplacian with Dirichlet boundary condition times a 
constant k > and thus its eigenpairs are given by 

ej(x) = y/2sin(jirx) and Xj = kir 2 j 2 for x G (0, 1), j G N. 

For each Fourier mode, we obtain from (I5.2p that 

- x h (Cm,e j ) H = V2e- x i h [ ( m (x)em(jirx)dx, j = l,2,...,N 

Jo 

(5.3) 



Y N ' M ' K t > 

1 m+l > c 3 I ~ ' 



18 



for m = 0, 1, M — 1 and 



(Y N ' M ' K , ej ) ={Z, ej ) H = V2 [ £(x)sm(j7rx)dx, j = l,2,...,N. (5.4) 
x ' H Jo 

Therefore, the implementation procedure goes as follows. Given Yj^' M,K , one can obtain 
Cm by (15. ip . Then we use some numerical integration method (here we choose com- 
posite trapezoidal formula) to approximate (( m ,ej) H in (15.31) for j = 1,2, ...,N. With 
(Cm,Gj) H ,j = 1,2, ...,N at hand, we can get Y^^f' K by (15.21) . Since the eigenfunc- 
tion 6j(x) = v / 2sin(j7rx) are sine functions, we can invoke built-in functions "idst" and 
"dst" in matlab to perform efficient computations. Recall that "dst" is a discrete sine 
transform, which transforms iV real numbers z(k),k = 1,2,...,N to A" real numbers 
y(j),j = 1, 2, N according to the following formula 

^ k 

fc=i 

and that the "idst" function is an inverse discrete sine transform, which transforms N 
real numbers z(k), k = 1, 2, N to N real numbers y(j),j = 1, 2, N according to the 
following formula 

2 " k 

y ti) = w^T[ z2 z ^ sin ^ 7r ArTT' ) ' 3 = lf '"' N ' 
k=i 

Setting z(k) = Cm(jv_b[) f° r ^ = 1,2, ...,iV in (15.61) . y(j) in ( 15. 6 p is in fact a composite 
trapezoidal formula to numerically approximate v2 (Cm, e i)- Hence 

(y^H' K ,e j ) H »e-^ h xy(j)/y/2, J = 1, 2, iV. (5.7) 

To be precise, only iV function values of Yj^' M ' K at grid points j^,k = 1,2, ...,iV are 
used to calculate A" function values of Cm at grid points jyTj, = 1, 2, iV. And then 
"idst" is used to approximate (Y m ]_ 1 ' , ej)n by (15.71) for j = 1,2,.. .,N. After that, we 
use a discrete sine transform "dst" to calculate A" function values of 5^4-1 ' at grid 
points, which are used to get A" function values of Cm at grid points before carrying 
out numerical integration at next step. Repeating this procedure we can finally obtain 
Y M ,M ' . In Figure [TJ we give the detailed implementation code of our new scheme ( II. 8p 
for the first test example. 

Before closing this subsection, we would like to give some remarks. Here and below, 
the aliasing errors caused by using composite trapezoidal formula are neglected and are 
not analyzed mathematically. As illustrated in the following numerical simulations, such 
errors do not effect the order of convergence. 



5.2 Numerical tests 

As the first numerical experiment, we consider an example SPDE ( II .ip in the introduction 
part with initial data £(x) = for x G (0, 1), T — 1, k — ^j, and 

f(x,y) = l- 2y, g{x,y) = - ^ 2 , fij = — , r)j(x) = e^x) = V2sm(jicx) (5.8) 
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N = 128; M = N~2; A = -pi~2*(l :N) . "2/200; mu=(l :N) . ~-2; 
f = @(x) l-2*x; g = @(x) (x+sin(x) ."3) ./((l+x.~2) .~2) ; 
Y = zeros (1 ,N) ; 

eta = zeros (1,N); SqrM=sqrt (M) ; 
for n=l:N 

eta = eta + 2*sin(n*(l :N)/(N+l)*pi) . ~2*mu(n) ; 

end 

for m = 1 :M 

y = dst(Y)*sqrt(2) ; 

dW = dst(randn(l,N) . *sqrt (mu*2/M) ) ; 

g_eva = g(y) ; 

y = y + f (y)/M+g_eva.*dW + . 5*SqrM*(g(y+g_eva/SqrM)-g_eva) . *(dW. ~2-eta/M) ; 
Y = exp(A/M) .*idst(y)/sqrt(2) ; 

end 

plot ( (0 : N+l) / (N+l) , [0 , dst (Y) *sqrt (2) , 0] ) ; 



Figure 1: Matlab code to simulate one path by the Runge-Kutta type scheme (11. 8p 
applied to the SPDE (11. ip with parameters as (15. 8p . 



for all x G (0, 1), y G R, j G N. Similarly to [HI Section 4], one can show that in this 
case the conditions (12~5) and (J277I) are fulfilled for 5 G (0, \),a G (0, §),7 G (|, |). As a 

•1 3 



result, Assumption 12. 1112751 are all satisfied for (3 = |,<5 G (0, G (0, |),7 G K > . i ,. 

According to the computational analysis in j!) , we know that the linear implicit Euler 
method (USD with M = N 3 ,K = N promises the existence of some real constants C r > 
and arbitrarily small r G (0, |) such that 



E 



X T -Y T 



N,N 3 ,N 



H 



(5.9) 



For the infinite version of Milstein type method (jl.7p . it is shown in [9] that, iV 2 time 
steps, in contrast to iV 3 time steps for the linear implicit Euler scheme ( II .5p . are required 
to achieve ( I5.9p . that is, (ll.7p with iV 2 time steps guarantees that for some real constants 
C r > and arbitrarily small r G (0, |) 



E 



Xt-Y, 



ryN,N 2 ,N 



N- 



H 



< C r N r 



(5.10) 



For the new scheme (11.81) applied to (II. ip . the main result (Theorem 12. ip in this paper 
shows that for some positive constants C r 



E 



Xt~Y m 



N.M.K 



2 \ 2 
H 



(5.11; 



holds for all arbitrarily small r G (0, |). Choosing M = N 2 ,K = N in the preceding 
result gives 



E 



for arbitrarily small < r < |. 



v V N,N 2 .N 
A T TV 2 



2 \ 2 



H 



< C r N V ~2 



(5.12) 
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In Figure [H we present detailed implementation code of the scheme ( 11.8J) . The term 
^ £7=1 ^jiVj) 2 i n (H -8p is computed once in advance for which 0(N 2 ) computational 
operations are needed. After that, 0(Nlog(N)) further computational operations and 
independent standard normal random variables are needed to compute ,N from 

Y^' n2 ' n by using the fast Fourier transform. Since iV 2 time steps are used, 0(N 3 log(iV)) 
computational operations and random variables are required to obtain Y^ N ' ■ Taking 
the convergence order |— in (I5.12p into account shows that the scheme (ll.8p promises the 
overall convergence order |— , which is the same as that of Milstein type scheme (jl.7p . 
But for the linear implicit Euler scheme f ll.5p . N 3 time steps are used and one can just 
get an overall convergence order of |— . 




Figure 2: Numerical results for SPDE (11.11) with parameters as (15. 8p . 

Figure [2] depicts approximation errors (II. 4p of the various approximations Y^ N ' N , 
Y$ N ' N ,Y^ N ' N with N = 4,8,16,32,64 against the number of used normal random 
variables on a log-log scale. As a measure for the computational effort, here we take the 
number of realizations of independent random variables needed for the calculation of the 
approximation. One can detect that the numerical results are consistent with our asser- 
tions on the convergence order. Besides, the Runge-Kutta method (11.81) and the Milstein 
type scheme (11.71) produce nearly the same approximation error. Numerical results also 
show that both the Runge-Kutta method (jl.8p and the Milstein type scheme (ll.7p are 
much more computationally effective than the linear implicit Euler scheme (ll.5p . For 
instance, Y^ 64 ' 6 in the case of the linear implicit Euler scheme and V^ 4 .' 64 ' 64 , Y^^f A ,64 
in the case of the Runge-Kutta method and the Milstein type scheme achieve the preci- 
sion e = 0.001 in (11.41) . For one path, it needs to generate 64 4 = 16777216 independent 
normal random variables and costs 105.312000 seconds to simulate Y^ m ' 64 . But the 
number of random variables needed to generate decreases to 64 3 = 262144 for one path 
simulation of Y^ M ' 64 , l^t' 64 ' 64 . Accordingly, the runtime for one path simulation of 

Y^ M ,64 and l^' 64 ,64 are, respectively, reduced to 2.328000 seconds and 1.984000 sec- 
onds. In Table (TJ we list runtime (seconds) for one path simulation using the three 
methods with various iV (N = 32,64,128,256). Note that the "exact" mild solution 
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is identified with the numerical solution using very small stepsize and that the matlab 
codes to simulate Y^ N ' N , Y^ N ,N and Y^ N ,N are presented in Figure 2, Figure 3 from 
[9] and Figure Q] in this article, respectively. It turns out that the schemes fll.8p and 
(11. 7p are progressively faster than the linear implicit Euler method (II. 5p as iV increases. 
Also, we observe from Table [T] that the Runge-Kutta type scheme (II. 8p is faster than 
the Milstein type scheme (II. 7p . This is due to the fact that evaluation of the partial 



derivative 
of the function g(x, y) 



3(l+i/ 2 )sin 2 (i/) cos(;/)-4i/siir i fa)-3'i/ 2 +l 



y+sm A (y) 

(i+y 2 ) 2 ' 



(i+y 2 ) 3 



costs more time than evaluation 



Table 1 


Runtime (seconds) for one 


path 


simulation using 


the three schemes 


yN,N 3 ,N 
1 m i 


Y£' N2 > N , Y£> n2 > n with N = 32, 64 


,128, 


256 






Linear Implicit Euler scheme 




Milstein scheme 


Runge-Kutta scheme 


N — 32 


10.172000 




0.547000 


0.469000 


N = 64 


105.312000 




2.328000 


1.984000 


N = 128 


1068.969000 




12.547000 


10.344000 


N = 256 


12604.844000 




77.016000 


60.312000 



As the second numerical experiment, we consider the case, where the two operators 
A and Q do not share the same eigenfunctions. More accurately, we choose the initial 
data £(x) = for x G (0, 1), T = 1, k — and the other parameters are set as 



f{x,y) = l-y,g(x,y) 



i + y 2 



Ho 



0, fj,j = — , r) (x) = 1, r)j(x) = V2cos(jirx) (5.13) 



for all x e (0,1), yeR, j G N. 

For this example, it is shown in [9] that Assumption 12.1112.51 are all satisfied with 
(3 = g,5 G (0, |),a G (0, |), 7 G Consequently, in this case Theorem 12.11 shows 

that 



E 



X T -Y, 



N,N 2 ,N 



iV 2 



H 



< C r N r - 



(5.14) 



holds for some positive constants C r , all arbitrarily small r G (0, 2) and iV G N. Hence 
its overall convergence order is |— . For the linear implicit Euler scheme (II. 5p . 



E 



X T (x) 



< C r N r ~ 



(5.15) 



holds for some positive constant C r , all arbitrarily small r G (0, 2) and iV G N. (I5.15P 
implies that the linear implicit Euler scheme has the overall convergence order of |— . 
These asymptotic results can be observed clearly in Figure EJ where approximation errors 
of the three approximations Y^ N ,N , Y^ N ,N , Y^ N ' N with = 2, 4, 8, 16, 32 against the 
number of used normal random variables are plotted. 
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